Precision benchmark calculations for four particles at unitarity 

Shahin Bour a , Xin U b , Dean Lee 6 , Ulf-G. Meifiner a ' c , Lubos Mitas fc 
a Helmholtz-Institut fur Strahlen- und Kernphysik (Theorie) and 
Bethe Center for Theoretical Physics, 
Universitat Bonn, D- 531 15 Bonn, Germany 
b Department of Physics, North Carolina State University, Raleigh, NC 27695, USA 
c Institut fur Kernphysik (IKP-3), Institute for Advanced Simulation (IAS-4), 
and Jiilich Center for Hadron Physics, 
Forschungszentrum Jiilich, D-52^25 Jiilich, Germany 

Abstract 

The unitarity limit describes interacting particles where the range of the interaction is zero and 
the scattering length is infinite. We present precision benchmark calculations for two-component 
fermions at unitarity using three different ab initio methods: Hamiltonian lattice formalism using 
iterated eigenvector methods, Euclidean lattice formalism with auxiliary- field projection Monte 
Carlo, and continuum diffusion Monte Carlo with fixed and released nodes. We have calculated 
the ground state energy of the unpolarized four-particle system in a periodic cube as a dimensionless 
fraction of the ground state energy for the non-interacting system. We obtain values 0.211(2) and 
0.210(2) using two different Hamiltonian lattice representations, 0.206(9) using Euclidean lattice, 
and an upper bound of 0.212(2) from fixed-node diffusion Monte Carlo. Released-node calculations 
starting from the fixed-node result yield a decrease of less than 0.002 over a propagation of OAEp 1 
in Euclidean time, where Ep is the Fermi energy. We find good agreement among all three ab 
initio methods. 
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I. INTRODUCTION 



The unitarity limit describes interacting particles where the range of the interaction is 
zero and the S-wave scattering length is infinite. In this paper we consider the unitarity 
limit of two-component fermions. Throughout our discussion we refer to the two degenerate 
components as up and down spins, though the correspondence with actual spin is not neces- 
sary. At sufficiently low temperatures the spin-unpolarized system is an S-wave superfluid 
with properties in between a Bardeen-Cooper-Schrieffer (BCS) fermionic superfluid at weak 
coupling and a Bose-Einstein condensate of dimers at strong coupling jl-3]. In nuclear 
physics the phenomenology of the unitarity limit approximately describes cold dilute neu- 
tron matter. The scattering length for elastic neutron-neutron collisions is about —18 fm 
while the range of the interaction is roughly the Compton wavelength of the pion, 1.4 fm. 
The unitarity limit is approximately realized when the interparticle spacing is about 5 fm. 
While these conditions cannot be produced experimentally, neutrons at around this density 
can be found in the inner crust of neutron stars. 

Experimental probes of the unitarity limit are now well established using trapped ul- 
tracold Fermi gases of alkali atoms. The characteristic length scale for the interatomic 
potential is the van der Waals length £ v dw- In the dilute limit the spacing between atoms 
can be made much larger than £ v dw an d the interatomic potential is well approximated 
by a zero-range interaction. The S-wave scattering length can be tuned using a magnetic 
Feshbach resonance 4 s|. This technique involves setting the energy level for a molecular 
bound state in a "closed" hyperfine channel to cross the scattering threshold for the "open" 
channel. The total magnetic moments for the two channels are different, and so the crossing 
can be produced using an applied magnetic field. 

The ground state for two-component fermions in the unitarity limit has no physical 
length scales other than the average distance between particles. The scaling properties in 
the unitarity limit are the same as that of a non-interacting Fermi gas. For Nf up spins 
and iVj, down spins in a given volume we write the energy of the unitarity-limit ground state 
as E^j N . For the same volume we call the energy of the free non- interacting ground state 
Eft ■ e £ ■ In the following we write the dimensionless ratio of the two energies as £n v n±, 

£ NtiNi = E^JE /*^. (1) 
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The parameter £ is defined as the thermodynamic limit for the spin-unpolarized system, 



£ = hm £ NN . 

iV->oo 



(2) 



II. RESULTS FOR £ AND THE NEED FOR PRECISION BENCHMARKS 

Several experiments have measured £ using the expansion rate of 6 Li and 40 K released 
from a harmonic trap as well as sound propagation. Some recent measured values for £ are 



0.32 
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<L 0.36(15) Q, 0.51(4) H, 0.46(5) [12], 0.46^ Q, 0.435(15) Q, 0.41(15 jH, 



0.41(2) 



16|, and 0.39(2) A new preliminary measurement finds a value 0.36(1) 



There are numerous analytical calculations of £ using a variety of techniques such as 



saddle point and variational approximations 18|, |l9) , Pade approximations and truncated 



series methods 



20J-I22J, mean field 



;rapolated from small systems 
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theory with pairing 



23 



241 ]. density functional theory ex- 



25j , renormalization group flow 



33], large- N expansions [34j, and other methods 35j. The values for £ range from 0.2 



26l | . dimensional expansions 



to 0.6 with most predictions in the range from 0.3 to 0.4. 

There are also many numerical calculations for £. The earliest fixed-node diffusion Monte 
Carlo simulations for N spin-up and N spin-down fermions in a periodic cube found £jv n to 



arger N 37J, [38] . A restricted path integral 



39] . and a sign- restricted mean field lattice 



be 0.44(1) for 5 < N < 21 36j and 0.42(1) for 
Monte Carlo calculation found similar results 

calculation yields 0.449(9) |40|. Another fixed-node diffusion Monte Carlo calculation sets 
an upper bound for £ N>N at 0.4244(1) for N = 33 and 0.4339(1) for iV = 64 [4lj. A more 



recent fixed-node calculation sets an upper bound for £jv jA r at 0.383(1) for iV between 2 and 
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42| | . This study includes an extrapolation to the zero-range limit and an analysis of shell 



effects using density functional theory. We note that methods such as fixed-node diffusion 
Monte Carlo provide only an upper bound for the ground state energy. An unbiased estimate 
for the ground state energy requires releasing the nodal constraint over a propagation time 
comparable to the diffusion time for neighboring particles to cross paths. 

There have also been a number of lattice simulations of two-component fermions in the 
unitarity limit. Several lattice simulations for the average energy at nonzero temperature 
have been extra pol ated to the zero temperature limit. The extrapolated zero temperature 
results from 3 established a bound, 0.07 < £ < 0.42. The results of Ref. 45] as well as 



Ref. 



461 . 1471 ] produce a value for £ in the 0.3 to 0.5 range. More recent lattice calculations 



extrapolated to zero temperature yield values £ = 0.292(24) [48|, |49| and £ = 0.37(5) [50 ]. 

In Ref. 5]J the ground state energy was calculated on the lattice using auxiliary-field 
Monte Carlo and Euclidean time projection starting from an initial state. The value of 
£iv,iv for iV = 3, 5, 7, 9, 11 were calculated at lattice volumes 4 3 , 5 3 , 6 3 in units of lattice 
spacing. From these small volumes it was estimated that £ = 0.25(3). In Ref. 52j this 
lattice calculation was improved using bounded continuous auxiliary fields. This calculation 
included an extrapolation to the continuum limit for £5 5 and £7.7 using lattice volumes 
4 3 , 5 3 , 6 3 , 7 3 , 8 3 . The results obtained were £5,5 = 0.292(12) and £7,7 = 0.329(5). Another 
technique called the symmetric heavy-light ansatz found similar values for £jv,at- While 
this approach is not an ab initio method, the agreement with the values for £55 and £77 in 



Ref. 



52] were within an error of 0.015. This method gives an estimate of £ = 0.31(1) in the 



continuum and thermodynamic limits [53]. Another extrapolation of the same data using 
density functional theory to include shell effects yields a value £ = 0.322(2) Some 
newer but preliminary lattice calculations using different projection and sampling methods 
produce a value £n,n — 0.412(4) for N in the range from 8 to 19 54], 1551 ] . 

The physics of the unitarity limit is universal and can be observed in many different 
systems and calculated using many different methods. However the spread in experimental, 
analytical, and numerical evaluations for £n,n and £ highlights the need for precision bench- 
marks and a more careful understanding of residual errors. Benchmarks at unitarity have 
been a subject of much discussion at several recent workshops and programs at the Institute 
for Nuclear Theory in Seattle. In this paper we discuss benchmarks for four unpolarized 
particles in a periodic cube. We focus on first principles numerical calculations for £2,2 
where all stochastic, extrapolation, and systematic errors can be reliably estimated. The 
three calculations we compare are the Hamiltonian lattice formalism using iterated eigenvec- 
tor methods, Euclidean lattice formalism with auxiliary-field projection Monte Carlo, and 
continuum diffusion Monte Carlo with fixed and released nodes. 



III. NOTATION AND DEFINITIONS 

Let Eft ™ be the ground state energy for up-spin and N± down-spin free fermions 
with equal masses in a periodic cube. We write E^ N for the ground state energy at 
unitarity for the same particle numbers, and N±, and the same periodic cube. In the 
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introduction we defined the energy ratio, 

C — P / E7i0,free fn\ 

We should point out that there are actually two different conventions for £jv tl jVj. used in the 
literature. We refer to Eq. d3J) as the few-body definition for the energy ratio 6^,7^- This 
is the definition we use for all calculations presented here. 

The alternative definition for the energy ratio £,N t ,N l is what we call the thermodynamical 
definition. This involves replacing E^™^ by the formula one gets in the thermodynamic 
limit. We define the Fermi momenta and energies in terms of the particle density, 

^=( 67r2 §) 1/3 ' ^=( 67r2 §) V3 ' ( 4 ) 

kp kp 
2m 2m 

In the thermodynamic limit the ground state energy of the non-interacting system is 

-N t E F j + -N^E^. (6) 

We use this to define the thermodynamical definition of the energy ratio, 

E° 

/■thermo _ N t> N l /y\ 

iN^ + lN.Ep,, {7) 

For finite and iVj. the few-body ratio £N t ,N± and thermodynamical ratio Cn^n° differ 
due to shell effects in the non-interacting system. There are several calculations in the 
literature using each of these two alternative definitions. In Table [J we have tabulated the 
conversion between the two definitions for several values of particle number with N-f = N±. 



IV. HAMILTON! AN LATTICE WITH SPARSE-MATRIX EIGENVECTOR IT- 
ERATION 



A. Formalism and notation 



Let n denote spatial lattice points on a three-dimensional L x L x L periodic cube. We 
use lattice units where physical quantities are multiplied by powers of the spatial lattice 
spacing to make the combination dimensionless. The two-component fermions are labelled 
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TABLE I: Conversion factor between the two ground state ratios Cat^a^ and £a^ c ™° for various 
values 7V| = N±. 
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1.0246 
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1.0064 



as spin- up and spin-down, the lattice annihilation operators are written as a^{n) and a±(n). 
We start with the free non-relativistic lattice Hamiltonian, 

2^ a\{n)Oi{n)- — 

n,j=t4 ' i=l,2,3 n,i=t4 

We define the spin-density operators 



i? free = — a}(n)tii(n) - ^ Ja}(n)oi(n + Z) + a\(n)ai(n - /)] . (8) 



= a|(n)a-f-(n), 



(9) 
(10) 



We consider two different lattice Hamiltonians each of which yield the unitarity limit in the 
low-energy limit. The first Hamiltonian, Hi, has a single-site contact interaction, 



Hi = H iTCC + Ciy]pt(n)p±(n). 



(11) 



The coefficient of C\ is tuned to set the S-wave scattering length a® to infinity. The second 
Hamiltonian, H 2 , has a contact interaction as well as nearest-neighbor interaction terms, 



H 2 = H hcc + C 2 ^ Pt(™)Pl(n) 



(12) 



J=l,2,3 n 

The coefficients C2 and C 2 are tuned so that ao goes to infinity while the S-wave effective 
range parameter r vanishes. 



6 



We use Liischer's formula 56M59| to determine the unknown interaction coefficients C\, 
C*2, and C' 2 . Liischer's formula relates the two-particle energy levels in a length L periodic 
cube to the S-wave phase shift, 

{ Lp \ 2 

V = 



1 



p cot 5 (p) = —S (77) , 

where S(rj) is the three-dimensional zeta function 

0(A 



Lp 
2^ 



lim 

A— >oo 



2 n 2 ) 



n z — 7] 



47rA 



In terms of 77, the energy of the two-particle scattering state is 



E- 



pole 



p 



m m 

The S-wave effective range expansion gives 



p cot 5 (p) = — - + \r Q p 2 + 0{p 4 ). 
a 2 



(13) 



(14) 



(15) 



(16) 



Setting a to infinity requires p cot S (p) to vanish at threshold. Setting both a to infinity 
and r to zero requires that pcot5o(p) is 0(p A ) near threshold. The plots for pcot5o(p) 
versus p 2 are shown in Fig. [TJ The values we find for the interaction coefficients are 



mCo 



-3.7235, 



-3.9570, 



mC 2 



-0.3008. 



(17) 
(18) 



B. Results for the four-particle benchmark 

n 

Using the Lanczos algorithm for sparse- matrix eigenvector iteration [60] . we have com- 
puted the ground state energy for two spin-up and two spin-down particles in a periodic 
cube of length L. For both lattice Hamiltonians, H\ and H 2 , we have computed £2,2 as 
defined in Eq. fl3J) for values L = 4,5, 6, 7, 8. The results are shown in Fig. (J2J). 

We have fitted the data using polynomials in 1/L up to third order and extrapolate to 
the infinite L limit with an estimated extrapolation error of ±0.002. We note that this 
extrapolation should remove all measurable lattice discretization effects. For Hi we find 

6, 2 = 0.211(2), (19) 
7 
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FIG. 1: (Color online) Plot of pcot<5o(p) versus p 2 for the lattice Hamiltonians Hi and Hi- 
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FIG. 2: (Color online) Ground state energy ratio £2,2 for lattice Hamiltonians H\ and i?2- We 
show results for values L = 4, 5, 6, 7, 8, and extrapolate to the infinite volume limit. 

and for H 2 we get 

6, 2 = 0.210(2). (20) 

The agreement between these two independent calculations is consistent with our estimate 
of the systematic errors. 

The third-degree polynomial extrapolation is made possible by the high precision data 
obtained for each L using Lanczos eigenvector iteration. For the Monte Carlo data appearing 
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later in our discussion we use only linear extrapolations in 1/L. For the H2 data we note 
the small slope in 1/L near 1/L = 0. This is expected due to the effective range ro being 
set to zero for H 2 - The small amount of linear dependence in 1/L that remains is likely 
due to other lattice artifacts such as the breaking of Galilean invariance 61]. 



V. EUCLIDEAN LATTICE WITH AUXILIARY-FIELD PROJECTION MONTE 
CARLO 

A. Formalism and notation 



For the Euclidean lat 
used in Ref. 



51 



Ace calculation we use the normal-ordered transfer-matrix formalism 



52 



62 



63] . Normal-ordering refers to the rearrangement of operators with 



annihilation operators on the right and creation operators on the left. This prescription 
is useful in that it provides an exact relation between Grassmann path integration and the 
operator formalism 



found in Ref. 



52 



62 



64 



65] . The details of the application of this correspondence can be 



63| . As before we use lattice units which correspond with multiplying 



physical quantities by the corresponding power of the spatial lattice spacing to make the 
combination dimensionless. We write at = a t /a for the dimensionless ratio of the temporal 
lattice spacing to spatial lattice spacing. For fermion mass m, we take the ratio of lattice 
spacings so that ra~ x a t = 0.1109. 

We start with the normal-ordered transfer matrix operator, 



M 



exp 



-H hee a t - Ca t ^ P\{n)Px{n) 



(21) 



The free lattice Hamiltonian Hf ree was defined in Eq. flH]) and the spin densities were defined 
Eq. f l9~fTUj) . Just as in the Hamiltonian lattice calculation we use Luscher's formula to 
determine the unknown coefficient C. Setting the S-wave scattering length to infinity, we 
find 

m C = -3.4938. (22) 

In Fig. ([3]) we plot pcot So(p) versus p 2 at unitarity. 

For the Monte Carlo simulations we use the bounded continuous auxiliary-field trans- 



formation introduced in Ref. 



52] . This transformation was shown in Ref. 52] to have 
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FIG. 3: (Color online) Plot of pcot5o(p) versus p 2 for the Euclidean-time lattice formalism. 



performance advantages over other continuous and discrete auxiliary-field transformations. 
We can write the transfer matrix as 



M = II ^r"/ ds (n,n t ) 



M(s,nt) 



where 



M(s, n t ) = : exp j-# frcc a t + ^ V~ 2Ca t sin [s(n, n t )} ■ [p t (rf) + p;(n)]| 



(23) 



(24) 



Let \^f l 22 ) a Slater-determinant of single-particle normal modes composed of two spin-up 
fermions and two spin-down fermions, 



|*2 n 2 t > = l^i)A|^ 2 )A|^3)A|^4). 
For the calculations presented here we choose 

= \/5 5 a t (T?) |0) ' ^ = \ZJ^ cos(2irn 3 /L)a\(n) |0> 



(25) 



(26) 



IV> 3 ) 



£ °l(n) |0> , |V>4> = / J E cos(27™ 3 /L)a|(n) |0> . (27) 



|V>i) and |-0 3 ) are constant-valued throughout the periodic box, while \ip 2 ) and l^) are 
standing waves with wavelength L along the z-axis. We construct the Euclidean-time 
projection amplitude 

" 1 



z,m s n [hj_ 

n,n t 



ds(n } nt) 



y^\M(s,L t -l) M(s,0)\¥$), (28) 
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TABLE II: Lattice dimensions L 3 x Lt used in calculations for £2,2 (t). 
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where the Euclidean projection time t equals L t times the temporal lattice spacing. 

Each M(s,n t ) consists of only single-particle operators interacting with the background 
auxiliary field. Therefore we find 



M(s, L t - 1) M(s,0) 1*$) = [detM(M)] 5 



(29) 



where 



[M(s, *)U = (Vv I L t - 1) M(s, 0) |^ fe ) , (30) 

for matrix indices k, k! = 1, 2, 3, 4. We define a t-dependent energy expectation value, 

(31) 



£! , 2 ( t ) = il„^('- a ') 



a* ^2,2 (t) 

We can also express the t-dependent expectation value as a fraction of the non-interacting 
ground state energy, 

^{t)=E 2 ^{t)/E^. (32) 
The ground state energy n ls gi ven by the limit 



El 2 = lim E 2 , 2 {t) 



t— >oo 



(33) 



and the desired few-body energy ratio can be computed as the limit 



£2,2 = lim 6,2 (t). 

1— >oo 



(34) 



B. Results for the four-particle benchmark 

For the calculation of £2,2^) we use the lattice dimensions L 3 x L t shown in Table HT1 
The simulations are run with 2048 processors each independently generating 3000 hybrid 
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FIG. 4: (Color online) The lattice data for £2,2(2) versus Ept for L = 4,5, 6, 7, 8. Also shown are 
the results of the asymptotic fits. 



Monte Carlo trajectories |66l-l68l|. To extract £2,2 we perform a least-squares fit of £2,2(2) to 
the asymptotic form, 

6, 2 (2)= 6,2 + be- SEt . (35) 
This exponential form takes into account the contribution from higher-energy states and 



is the same method used in Ref. 



51 



52|, 



62 



63| . We focus on measuring the ground 



state energy accurately and ignore the numerically small contributions hidden in the far 
asymptotic tail of £2,2(2). We determine b, 5E, and £2,2 from least-squares fitting over the 
range Ept = 2 to Ept = 6. The lattice data for £2,2(2) together with the asymptotic fits are 
shown in Fig. HI 

We use the lattice results for £ 2> 2 with L = 4,5, 6, 7, 8 to extrapolate to the continuum 
limit L — > 00. In Fig. |5]we show the lattice results for £ 2 ^ for L — 4,5, 6, 7, 8 plotted versus 
L~ l . We expect a dependence on L arising from effects such as the effective range correction 
and lattice cutoff effects. Using a linear extrapolation in L" 1 , we obtain the continuum 
limit value 

£2,2 = 0.206(9). (36) 

This result using Euclidean lattice projection Monte Carlo is in agreement with the Hamil- 
tonian lattice results in Eq. (I19|20p . 
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FIG. 5: (Color online) Results for £2,2 for L = 4,5,6,7,8 plotted versus L _1 . The lattice results 
are extrapolated to the continuum limit L — > 00. 

VI. DIFFUSION MONTE CARLO WITH FIXED AND RELEASED NODES 
A. Formalism and notation 

We now discuss diffusion Monte Carlo calculations for the same benchmark system of 
four particles at unitarity in a periodic cube with length L. This time, however, we use 
continuous variables and consider the evolution as a function of Euclidean time as a diffusion 
Monte Carlo process using an ensemble of random walkers. An introduction to the basic 
techniques can be found in Ref. 69_|. For the interaction between spin-up and spin-down 
fermions we use a Poschl- Teller potential tuned to infinite S-wave scattering length. For 
fermion mass m, the form of the potential is 



V{r) 



m cosh 2 (^r) ' 

where the momentum scale /i determines the S-wave effective range parameter, 

r = 2fi~\ 



(37) 



(38) 



The unitarity limit corresponds with taking the limit fi — > oo. For the unpolarized four- 
particle system we let R be the set of individual particle coordinates, 



R= {n v r lv r 2v r 2l } 



(39) 
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We use a BCS-type pairing wavefunction projected onto two spin-up and two spin-down 
fermions. The wavefunction ^bcs(R') can be written as a 2 x 2 Slater determinant, 



BCS 



det 



(40) 



(r 2t ~ r h ) 0(f 2t - f 2 J 

where </>(r) is the pairing function. The trial wavefunction, \l/r(R), is given as a product of 
^bcs(R-) times a Jastrow factor |70[], 



* T (R) = *Bcs(R)exp[J(R)]. 



(41) 



The Jastrow factor incorporates particle correlations and is useful in reducing stochastic 
errors in the Monte Carlo calculation 7lM74| . We use a product of Gaussian functions for 



the Jastrow factor. The exponents of these Gaussian functions are tuned to minimize the 
combination of the variational energy and the variance of the local energy in the variational 
Monte Carlo. We note that the positive definite function exp [J(R)] has no effect on the 
nodal structure of Sl/WR). 



To determine the pairing function 0(r) we use an approach similar to that in Ref. 36]. 
We use an ansatz which is a superposition of Gaussian functions with periodic copies, 

0(f) = ^4 e^ +s ^ 2 e-^y +s ^ 2 e-^ z+s ^ 2 . (42) 



Here r = (x,y,z), and dk and at are variational parameters. Gaussian functions from 
the nearest periodic images at distance L make small contributions while periodic copies 
further away can be neglected. This construction has the advantage of providing more 
flexibility for pairing orbitals in the box while keeping the orbitals smooth over the periodic 
boundary with zero derivative. The Jastrow factor exp[J(R)] is constructed similarly. The 
calculation is performed using standard variational Monte Carlo sampling of the square of 
the trial function 



69 



75|. 



For the diffusion Monte Carlo calculation we use a large ensemble of forward-propagating 
random walkers with population branching processes determined by the local energy and 
guided by the trial wavefunction \&r(R). The local density of random walkers gives a 
statistical estimate of the product of the propagated quantum wavefunction \I/(R) times the 
\1/t(R) trial function. In the fixed-node diffusion Monte Carlo (FN-DMC) calculation, the 
nodal structure of ^(R) is fixed by ^(R), and the product \I/(R) ^(R) is therefore positive 
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definite. Since the trial function is known explicitly and analytically it is then possible to 
extract the energetics of the lowest fermionic state within the fixed-node boundary conditions 
69, t3- 

The fixed-node calculation sets an upper bound on the ground state energy To measure 
the quality of the upper bound we release the nodal constraints 76). For the released- node 
diffusion Monte Carlo calculation (RN-DMC) we use a positive-definite guiding profile for the 
diffusion of random walkers. In the calculations presented here we consider a one-parameter 
family of guiding profiles, 



*S(R) = v /^(R) + «(^). (43) 

The dimensionless parameter a controls the rate of diffusion across the nodal boundaries of 
^t(R-); an d (&t) i s the average value of \E^(R ) evaluated over all R , where R is the 
configuration right after the nodal release process. 



B. Results for the four-particle benchmark 

In Fig. [6] we show fixed-node (FN-DMC) and released-node (RN-DMC) diffusion Monte 
Carlo results for two spin-up and spin-down fermions in a periodic cube at unitarity. The 
ratio of the effective range parameter to the length of the cube is r^/L = 1/40. We use 
parameters a = 0.01, 0.05, 0.2 for the guiding profile \I/^(R). With each nodal crossing 
the associated weights of the random walkers change sign, leading to a sign cancellation 
problem which grows exponentially with Euclidean propagation time t. For the calculations 
presented here we have measured the released-node correction starting from the fixed-node 
result up to propagation time of t = OAEp . As seen in Fig. [6J the decrease in £2,2 is less 
than 0.002 over the duration of the released-node time propagation. 

We have repeated the fixed-node and released-node calculations of £2,2 for values r /L = 
1/20, 1/40, 1/80, 1/160, 1/320, 1/640. The fixed-node results for £ 2 ,2 versus the r /L are 
shown in Fig. [7J Using a linear fit in r Q /L, we extrapolated to the limit of zero effective 
range. 

In the zero-range limit we get the final result 

6,2 = 0.212(2) (44) 
in the fixed-node approximation. The error estimate includes the stochastic error, zero- 
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FIG. 6: (Color online) Fixed-node and released-node diffusion Monte Carlo results for tq/L = 0.025 
and a = 0.01, 0.05, 0.2. 
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FIG. 7: (Color online) Fixed-node diffusion Monte Carlo results for £2,2 versus the vq/L. The data 
is extrapolated to the limit of zero effective range. 

range extrapolation error, and errors due to time step discretization and other small residual 
effects. The released-node calculations are quantitatively similar to the results shown in 
Fig. [6] for tq/L = 1/40. The decrease in £ 2 ,2 is less than 0.002 over a propagation time 
of OAEp . We note that Ep 1 is the characteristic time scale required for two neighboring 
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particles of the same spin to cross paths. Given the stochastic noise in the released-node 
calculation results, it is difficult to pin down a correction to the fixed-node result from the 
nodal release. However a reasonable conservative estimate is that the upper bound set by the 
fixed-node calculation is less than 0.002/0.4 = 0.005 above the actual value. This appears 
confirmed by the agreement of Eq. (j44j) with the values for 6,2 obtain using Hamiltonian 
lattice eigenvector iteration and Euclidean lattice Monte Carlo. 



VII. SUMMARY AND DISCUSSION 



We have presented benchmark calculations for four unpolarized particles in a periodic 
cube using three different methods. In the Hamiltonian lattice formalism with iterated 
eigenvector methods, we obtained 

6,2 = 0.211(2) (45) 
using the Hamiltonian H\ defined in Eq. pip , and 

6,2 = 0.210(2) (46) 

using the Hamiltonian H 2 defined in Eq. (1121) . Using the Euclidean lattice formalism with 
auxiliary-field projection Monte Carlo we found the result 

6,2 = 0.206(9). (47) 

With fixed-node diffusion Monte Carlo in the continuum we extracted the upper bound 

6,2 < 0.212(2). (48) 

The release-node Monte Carlo calculation shows a decrease in 6,2 that is less than 0.002 
over a propagation time of OAEp 1 . We estimate that the upper bound set by the fixed- 
node calculation is less than 0.005 above the actual value. All three methods agree within 
estimated errors. The unpolarized four-particle benchmarks presented here should be useful 
for testing and calibrating residual errors for other numerical methods and perhaps also 
analytical calculations. We note that the comparison requires using the few-body definition 
of 6,2 in Eq. fl3]). If the thermodynamical definition in Eq. (JTj) is used then the conversion 
factor is 6,2 = 0.7331^ h 2 ermo . 
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We note the importance of continuum limit extrapolations for lattice calculations and 
zero-range limit extrapolations for continuum diffusion Monte Carlo calculations. The 
importance of continuum limit extrapolations in lattice calculations was already noted and 



measured in Ref. 



521 ] . In the fixed- node diffusion Monte Carlo calculations presented here 



we find that r /L needs to be less than 0.03 in order to obtain a value for £ 2 ,2 accurate to 
a relative error of 0.1. This is consistent with the range dependence found in Ref. 421. 

n 

After discussion with the authors of Ref. [42| , we are informed that they obtained an upper 
bound for £2,2 from fixed-node diffusion Monte Carlo agreeing within two significant digits 
with the results reported here. 

In this paper the four-particle benchmark was chosen to allow comparisons among several 
very different numerical methods. For larger N — — N± systems it is unfortunately not 
possible to use iterated eigenvector methods due to exponential L 6N ~ 3 scaling in memory. 
However the Euclidean lattice Monte Carlo and diffusion Monte Carlo methods extend 
readily to larger systems, and it is useful to comment on the relationship between the four- 
particle calculations discussed here and calculations in larger systems at unitarity. 

The four-particle results presented here should be useful for comparisons of different lat- 
tice Monte Carlo calculations using different lattice actions and algorithms. We note that 
there is a significant correction produced by the extrapolation to the continuum limit. For 
some lattice actions this extrapolation decreases the ground state energy while for others it 
increases the ground state energy. In all cases it is important that stochastic and systematic 
errors are sufficiently small for each chosen lattice spacing so that the continuum extrapo- 
lation can be done accurately. The same lattice methods used here to find £2.2 = 0.206(9) 
were also used to determine £ 5 5 = 0.292(12) and £7,7 = 0.329(5) in Ref. [52j. In order 
to benchmark different lattice Monte Carlo methods, one starting point would be to test 
agreement with each of these values for £n,n- 

After completion of the original manuscript of this paper, the fixed-node diffusion Monte 
Carlo methods presented here were also applied to larger N = iV-f = N± systems at unitarity. 
In Ref. H the results £7,7 < 0.407(2), £ 19i i 9 < 0.409(3), and £33,33 < 0.398(3) are presented. 



These results are comparable to the upper bounds found in Ref. [42J, and in each case a 
significant reduction in the ground state energy is seen when extrapolating to the zero range 
limit. For the released-node calculations in these larger systems, the exponential severity 
of sign cancellations makes it difficult to extract data for Euclidean propagation time t 
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greater than N~ 1 Ep 1 . As noted above, Ep 1 is the characteristic time scale required for 
two neighboring particles of the same spin to cross paths. For Af = 2 we extrapolated 
to get a bound on the decrease in £2,2 over a propagation time of Ep 1 . For larger N the 



We note that there remains a significant gap between the lattice result £^7 = 0.329(5) and 
the fixed-node upper bound £7,7 < 0.407(2). This discrepancy must be better understood. 
A good starting point would be to make benchmark comparisons for other small systems, 
N = 3,4, 5. First it should be established that different lattice calculations agree on the 
values for £ 3)3 , £ 4)4 , and £5,5. Next the difference between lattice and fixed-node results 
should be measured for each N. The key question then is if this difference can be resolved 
by released- node calculations for smaller values of N. Given the agreement for £ 2 ,2 ; there 
is some reason to suggest that this may be possible. 
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